Main Results
We now look at the main results of the experiment. We are going to look first at the effect of richness, temperature and nutrients on community temporal stability. Then, we are going to look at the effect of fundamental and realised balance on temporal stability. Finally, we are going to look at the relationship between response diversity and temporal stability.
In the whole analysis, we calculated the temporal stability of total community biomass as the inverse of the coefficient of variation (ICV) (i.e. \(\frac{\sigma}{\mu}\)).
Effect of T, N and R
We can see that richness does not have a clear effect on community temporal stability, while stability was higher at lower temperature, and nutrients increased community temporal stability.
Effect of Divergence
We look at the relationship between divergence (our original response diversity metric)
Divergence is negatively related to temporal stability, suggesting that response diversity promotes stability. However, the relationship between divergence and stability becomes weaker as richness increases. This is why, after running the experiment, we developed another metric to measure response diversity, which we called balance, and that is presented in the main text of the publication.
Here, we provide extensive evidence of why balance is a better metric to measure response diversity than divergence, and thus justifying focusing the analysis around balance.
Comparing Divergence and Balance
We fist compare the predictive power of divergence vs balance.
Balance
#
mod1 <- lm(data=complete_aggr,log10(stability)~log10(balance_f))
check_model(mod1)
performance(mod1)
## # Indices of model performance
##
## AIC | AICc | BIC | R2 | R2 (adj.) | RMSE | Sigma
## ---------------------------------------------------------------
## -89.273 | -89.173 | -78.794 | 0.192 | 0.188 | 0.151 | 0.152
Divergence
mod2 <- lm(data=complete_aggr,log10(stability)~(divergence))
check_model(mod2)
performance(mod2)
## # Indices of model performance
##
## AIC | AICc | BIC | R2 | R2 (adj.) | RMSE | Sigma
## ---------------------------------------------------------------
## -55.716 | -55.615 | -45.237 | 0.072 | 0.068 | 0.162 | 0.163
Divergence and balance are both negatively related to stability, but balance explains more of the variance in stability than divergence.
Moreover, from the plot above, it looks like divergence declines in performance as richness increases. Let’s test this analyticaly.
# getting model estimates for each richness level
lm_divergence_richness_E <- complete_aggr %>%
nest(data = -richness) %>%
mutate(
model = map(data, ~ lm(log10(stability) ~ scale(divergence), data = .x)),
results = map(model, broom::tidy)
) %>%
unnest(results) %>% dplyr::filter(term=="scale(divergence)")
# getting model R squared for each richness level
lm_divergence_richness_R <- complete_aggr %>%
nest(data = -richness) %>%
mutate(
model = map(data, ~ lm(log10(stability) ~ scale(divergence), data = .x)),
results = map(model, broom::glance)
) %>%
unnest(results)
lm_divergence_richness_R
## # A tibble: 3 × 15
## richness data model r.squared adj.r.squared sigma statistic p.value df
## <fct> <list> <list> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2 <tibble> <lm> 0.201 0.191 0.209 19.9 2.66e-5 1
## 2 3 <tibble> <lm> 0.172 0.161 0.108 16.4 1.19e-4 1
## 3 4 <tibble> <lm> 0.0337 0.0215 0.129 2.76 1.01e-1 1
## # ℹ 6 more variables: logLik <dbl>, AIC <dbl>, BIC <dbl>, deviance <dbl>,
## # df.residual <int>, nobs <int>
# getting model estimatesf or each richness level
lm_balance_richness_E <- complete_aggr %>%
nest(data = -richness) %>%
mutate(
model = map(data, ~ lm(log10(stability) ~ scale(log10(balance_f)), data = .x)),
results = map(model, broom::tidy)
) %>%
unnest(results) %>% dplyr::filter(term=="scale(log10(balance_f))")
summary(lm(log10(stability) ~ scale((divergence)), data = complete_aggr))
##
## Call:
## lm(formula = log10(stability) ~ scale((divergence)), data = complete_aggr)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.41454 -0.08948 -0.00742 0.09164 0.65272
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.24267 0.01042 -23.278 < 2e-16 ***
## scale((divergence)) 0.04520 0.01045 4.327 2.22e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1625 on 241 degrees of freedom
## Multiple R-squared: 0.07208, Adjusted R-squared: 0.06823
## F-statistic: 18.72 on 1 and 241 DF, p-value: 2.218e-05
# getting model R squared for each richness level
lm_balance_richness_R <- complete_aggr %>%
nest(data = -richness) %>%
mutate(
model = map(data, ~ lm(log10(stability) ~ scale(log10(balance_f)), data = .x)),
results = map(model, broom::glance)
) %>%
unnest(results)
lm_balance_richness_R
## # A tibble: 3 × 15
## richness data model r.squared adj.r.squared sigma statistic p.value df
## <fct> <list> <list> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2 <tibble> <lm> 0.188 0.177 0.210 18.2 5.36e-5 1
## 2 3 <tibble> <lm> 0.232 0.222 0.104 23.8 5.44e-6 1
## 3 4 <tibble> <lm> 0.272 0.263 0.112 29.6 5.84e-7 1
## # ℹ 6 more variables: logLik <dbl>, AIC <dbl>, BIC <dbl>, deviance <dbl>,
## # df.residual <int>, nobs <int>
We can see that the relationship between divergence and stability becomes weaker as richness increases, while the relationship between balance and stability remains stable across richness levels.
Now we build a linear model were stability is modelled as a function of balance and divergence.
Then, we compared the variance explained by the full model compared to a model containing either only balance or only divergence.
lm_div_balance <- lm(data=complete_aggr,log10(stability)~log10(balance_f)+divergence)
check_model(lm_div_balance)
performance(lm_div_balance)
## # Indices of model performance
##
## AIC | AICc | BIC | R2 | R2 (adj.) | RMSE | Sigma
## ---------------------------------------------------------------
## -88.977 | -88.809 | -75.005 | 0.197 | 0.191 | 0.151 | 0.151
model with only divergence
lm_div <- lm(data=complete_aggr,log10(stability)~divergence)
check_model(lm_div)
performance(lm_div)
## # Indices of model performance
##
## AIC | AICc | BIC | R2 | R2 (adj.) | RMSE | Sigma
## ---------------------------------------------------------------
## -55.716 | -55.615 | -45.237 | 0.072 | 0.068 | 0.162 | 0.163
model with only balance
lm_balance <- lm(data=complete_aggr,log10(stability)~log10(balance_f))
check_model(lm_balance)
## # Indices of model performance
##
## AIC | AICc | BIC | R2 | R2 (adj.) | RMSE | Sigma
## ---------------------------------------------------------------
## -89.273 | -89.173 | -78.794 | 0.192 | 0.188 | 0.151 | 0.152
anova(lm_div_balance, lm_balance)
## Analysis of Variance Table
##
## Model 1: log10(stability) ~ log10(balance_f) + divergence
## Model 2: log10(stability) ~ log10(balance_f)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 240 5.5044
## 2 241 5.5432 -1 -0.038724 1.6884 0.1951
anova(lm_div_balance, lm_div)
## Analysis of Variance Table
##
## Model 1: log10(stability) ~ log10(balance_f) + divergence
## Model 2: log10(stability) ~ divergence
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 240 5.5044
## 2 241 6.3640 -1 -0.85959 37.479 3.741e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Balance explains more of the variance in stability than divergence, and there is virtually no difference between a model containing only balance and the full model.
Richness had to transformed to numeric and to be centered to avoid collinearity with divergence
lm_rich_div <- lm(data=complete_aggr,log10(stability)~divergence*scale(as.numeric(richness)))
check_model(lm_rich_div)
anova(lm_rich_div)
## Analysis of Variance Table
##
## Response: log10(stability)
## Df Sum Sq Mean Sq F value Pr(>F)
## divergence 1 0.4943 0.49435 19.8282 1.301e-05 ***
## scale(as.numeric(richness)) 1 0.1558 0.15579 6.2487 0.013100 *
## divergence:scale(as.numeric(richness)) 1 0.2496 0.24958 10.0106 0.001758 **
## Residuals 239 5.9587 0.02493
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Divergence significantly interact with richness, suggesting that the relationship between divergence and stability changes with richness.
While an ideal metric of response diversity should be independent of richness.
We repeat the same model using balance instead of divergence.
lm_rich_balance <- lm(data=complete_aggr,log10(stability)~log10(balance_f)*scale(as.numeric(richness)))
check_model(lm_rich_balance)
anova(lm_rich_balance)
## Analysis of Variance Table
##
## Response: log10(stability)
## Df Sum Sq Mean Sq F value
## log10(balance_f) 1 1.3152 1.31522 57.1274
## scale(as.numeric(richness)) 1 0.0003 0.00028 0.0122
## log10(balance_f):scale(as.numeric(richness)) 1 0.0405 0.04050 1.7589
## Residuals 239 5.5024 0.02302
## Pr(>F)
## log10(balance_f) 8.694e-13 ***
## scale(as.numeric(richness)) 0.9123
## log10(balance_f):scale(as.numeric(richness)) 0.1860
## Residuals
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Balance does not significantly interact with richness, suggesting that the relationship between balance and stability is stable across richness levels.
Finally, we assess variable importance using the relative importance of predictors in the full model.
We use the package vip (https://cran.r-project.org/web/packages/vip/vignettes/vip.html) to calculate the relative importance of predictors in the full model.
The function vip::vip for multiple linear regression, or linear models (LMs), uses the absolute value of the 𝑡
-statistic as a measure of VI.7. Motivation for the use of the associated 𝑡-statistic is given in Bring (1994) [https://www.tandfonline.com/doi/abs/10.1080/00031305.1994.10476059].
vip::vip(lm_div_balance)
We believe that the extensive evidence here provided justifies focusing the analysis around balance, and not divergence, as a metric of response diversity.
We will thus only look at balance for the rest of the analysis.
Effect RD
We are now going to look at how response diversity (balance) affected temporal stability of total community biomass. We are going to look at the relationship between fundamental baalance (so based only on species response surfaces measured in monoculture), an realised balance (measured accounting for species contribution to balance).
This is fundamentally testing our most important hypothesis.
We can see that balance is always negatively related to temporal stability, which means that response diversity promotes stability across richness levels. Interestingly, we see that there is little difference between fundamental and realised balance. Yet, as the richness increases, the relationship between realised balance and stability becomes steeper compared to fundamental balance.
We can see that the positive relationship between temporal stability and response diversity measured as divergence holds, but it becomes shallower as richness increases. We could speculated that this due to divergence considering only the responses of the 2 most “responding” species. Thus, when species richness increases, disregarding the responses of the other species in the community except the 2 responding the most makes the relationship between response diversity and stability weaker.
Linear models
Model: Fundamental balance
First we analyze the effect of fundamental balance, temperature, nutrients and richness on biomass temporal stability using a linear model.
balance and richness were modelled as continuous variables, while temperature and nutrients were modelled as categorical variables. balance was log-transformed to meet the assumptions of linear models.
lm_full<-lm(data=complete_aggr,log10(stability)~log10(balance_f)+as.numeric(richness)+nutrients+temperature)
Check model’s assumptions
##
## Call:
## lm(formula = log10(stability) ~ log10(balance_f) + as.numeric(richness) +
## nutrients + temperature, data = complete_aggr)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.28605 -0.07635 -0.01196 0.04542 0.42770
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.350905 0.033407 -10.504 < 2e-16 ***
## log10(balance_f) -0.050863 0.016041 -3.171 0.00172 **
## as.numeric(richness) -0.006457 0.009431 -0.685 0.49425
## nutrients0.35 g/L 0.179918 0.018750 9.596 < 2e-16 ***
## nutrients0.75 g/L 0.213112 0.019473 10.944 < 2e-16 ***
## temperature22-25 °C -0.078022 0.018683 -4.176 4.18e-05 ***
## temperature25-28 °C -0.100991 0.024737 -4.083 6.10e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1186 on 236 degrees of freedom
## Multiple R-squared: 0.5158, Adjusted R-squared: 0.5035
## F-statistic: 41.9 on 6 and 236 DF, p-value: < 2.2e-16
## Analysis of Variance Table
##
## Response: log10(stability)
## Df Sum Sq Mean Sq F value Pr(>F)
## log10(balance_f) 1 1.3152 1.31522 93.4656 < 2.2e-16 ***
## as.numeric(richness) 1 0.0003 0.00028 0.0199 0.8879
## nutrients 2 1.8841 0.94203 66.9447 < 2.2e-16 ***
## temperature 2 0.3379 0.16896 12.0071 1.081e-05 ***
## Residuals 236 3.3209 0.01407
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A linear model was fitted to examine the effects of resource balance, richness, nutrients, and temperature on community stability (measured as log₁₀(stability)). The model explained a significant portion of the variance (Adjusted R² = 0.5115, F(7, 235) = 37.2, p < 2.2e-16).
The intercept of the model was estimated at -0.349 (SE = 0.028, p < 2e-16), indicating the baseline log₁₀(stability) when all predictor variables are at their reference levels.
Among the predictors, log₁₀(balance) showed a significant negative effect on stability (Estimate = -0.054, SE = 0.016, p = 0.0009). This suggests that as balance increases (more balance), stability tends to decrease.
Nutrient concentration also had a significant positive effect on stability, with estimates for 0.35 g/L (Estimate = 0.180, SE = 0.019, p < 2e-16) and 0.75 g/L (Estimate = 0.212, SE = 0.019, p < 2e-16) indicating increased stability with higher nutrient levels.
Finally, temperature regimes showed a significant effect on stability. Both 22–25 °C (Estimate = -0.078, SE = 0.019, p = 3.81e-05) and 25–28 °C (Estimate = -0.098, SE = 0.025, p = 8.44e-05) significantly reduced stability when compared to the baseline (18–21 °C).
Richness did not show a significant effect on stability (Estimate = 0.002, SE = 0.019, p = 0.91), which suggests that it is not richness per se that affects stability in this system. Similar results have been previously reported in the literature, e.g. (Petchey et al. 2002)[https://nsojournals.onlinelibrary.wiley.com/doi/full/10.1034/j.1600-0706.2002.990203.x?casa_token=THaSxpjziQcAAAAA%3Ay_0gJhnL_rcPsrolHEmZvI0VF14a43WRDTy_UB_kDPQOxq2EA98NQT65Co63J58NCeW-SiTIjoDkgelQ] and (Gonzalez and Descamps-Julien 2004) [https://nsojournals.onlinelibrary.wiley.com/doi/full/10.1111/j.0030-1299.2004.12925.x?casa_token=4OtlkUcK3YgAAAAA%3ASFjll62wMQW8CBV9qDe5lbGGDXVcRiWna6slN1GaVvTMSRJKDAnpSd_jEBsFKKRsL7muxbJ0Tv3T] who found that species richness did not significantly impact temporal stability in fluctuating environment.
In summary, our findings show that temporal stability is significantly influenced by response diversity (balance), nutrient concentration, and temperature, with higher nutrient concentrations enhancing stability and higher temperatures reducing it. However, species richness was not a significant determinant of stability within the conditions of this study.
Prepare publication-ready table
Summary table
| Characteristic |
Beta |
95% CI |
p-value |
| log10(balance_f) |
-0.05 |
-0.08, -0.02 |
0.002 |
| as.numeric(richness) |
-0.01 |
-0.03, 0.01 |
0.5 |
| nutrients |
|
|
|
| 0.01 g/L |
— |
— |
|
| 0.35 g/L |
0.18 |
0.14, 0.22 |
<0.001 |
| 0.75 g/L |
0.21 |
0.17, 0.25 |
<0.001 |
| temperature |
|
|
|
| 18-21 °C |
— |
— |
|
| 22-25 °C |
-0.08 |
-0.11, -0.04 |
<0.001 |
| 25-28 °C |
-0.10 |
-0.15, -0.05 |
<0.001 |
sum vs. weighted_sum
We now look at how fundamental and realised balance are related to each other.
SEM
Now, we use a structural equation model (SEM) to explore how stability is influenced by asynchrony, temperature, nutrient levels, balance, and richness, with asynchrony also modeled as dependent on balance, nutrients, and richness.
## lavaan 0.6-19 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 12
##
## Number of observations 241
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 1.777 1.537
## Degrees of freedom 1 1
## P-value (Chi-square) 0.183 0.215
## Scaling correction factor 1.156
## Satorra-Bentler correction
##
## Model Test Baseline Model:
##
## Test statistic 295.357 353.081
## Degrees of freedom 9 9
## P-value 0.000 0.000
## Scaling correction factor 0.837
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.997 0.998
## Tucker-Lewis Index (TLI) 0.976 0.986
##
## Robust Comparative Fit Index (CFI) 0.998
## Robust Tucker-Lewis Index (TLI) 0.981
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) 152.706 152.706
## Loglikelihood unrestricted model (H1) NA NA
##
## Akaike (AIC) -281.411 -281.411
## Bayesian (BIC) -239.594 -239.594
## Sample-size adjusted Bayesian (SABIC) -277.631 -277.631
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.057 0.047
## 90 Percent confidence interval - lower 0.000 0.000
## 90 Percent confidence interval - upper 0.192 0.176
## P-value H_0: RMSEA <= 0.050 0.306 0.356
## P-value H_0: RMSEA >= 0.080 0.531 0.469
##
## Robust RMSEA 0.051
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.200
## P-value H_0: Robust RMSEA <= 0.050 0.327
## P-value H_0: Robust RMSEA >= 0.080 0.525
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.012 0.012
##
## Parameter Estimates:
##
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## stability ~
## asynchrny_Grss 0.195 0.033 5.947 0.000 0.195 0.399
## temperature -0.049 0.011 -4.454 0.000 -0.049 -0.241
## nutrients 0.147 0.011 13.216 0.000 0.147 0.719
## log_balance_f -0.029 0.012 -2.452 0.014 -0.029 -0.120
## richness 0.011 0.010 1.181 0.238 0.011 0.055
## asynchrony_Gross ~
## log_balance_f -0.080 0.029 -2.719 0.007 -0.080 -0.159
## nutrients -0.210 0.022 -9.373 0.000 -0.210 -0.501
## richness -0.102 0.023 -4.404 0.000 -0.102 -0.243
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .stability -0.407 0.040 -10.122 0.000 -0.407 -2.441
## .asynchrny_Grss 0.125 0.065 1.921 0.055 0.125 0.367
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .stability 0.012 0.001 10.670 0.000 0.012 0.430
## .asynchrny_Grss 0.081 0.009 8.917 0.000 0.081 0.694
##
## R-Square:
## Estimate
## stability 0.570
## asynchrny_Grss 0.306
Model Fit
The model fit indices suggest a good fit:
Comparative Fit Index (CFI) = 0.998 and Tucker-Lewis Index (TLI) = 0.986, both indicating a good fit as values close to 1 are considered strong.
Root Mean Square Error of Approximation (RMSEA) = 0.047 (with robust RMSEA at 0.051) and Standardized Root Mean Square Residual (SRMR) = 0.012. These values indicate a good fit, with RMSEA and SRMR values below 0.05 generally preferred.
Interpretation of Pathways
Stability:
Asynchrony: Positive and highly significant effect on stability, suggesting that asynchrony (indicating lack of synchrony or compensatory dynamics) is associated with greater stability.
Temperature: Negative effect, where higher temperature values correlate with lower stability, potentially due to physiological stress or disruption in community dynamics at higher temperatures.
Nutrients: Positive and highly significant, suggesting that greater nutrient availability enhances stability, possibly through support for higher productivity or resource buffering.
balance: Negative and significant effect, where greater balance reduces stability.
Richness: Not significant, indicating that within this model, richness does not have a notable effect on stability.
Asynchrony:
balance: Negative and significant, suggesting that greater balance reduces asynchrony.
Nutrients: Negative and highly significant effect, indicating that higher nutrient concentrations are associated with lower asynchrony, possibly due to homogenizing effects of nutrient availability.
Richness: Negative and significant, where increased richness is associated with reduced asynchrony, possibly indicating increased interactions or overlap in resource use among species.
Explained Variance
Stability: The model explains 57% of the variance in stability, suggesting a substantial amount of stability is accounted for by these factors.
Asynchrony: The model explains 30.6% of the variance in asynchrony, indicating that while balance, nutrients, and richness contribute, other factors may also play a role in driving asynchrony.
Summary
This SEM model demonstrates that stability in the ecosystem is positively associated with asynchrony and nutrient levels, but negatively associated with temperature and balance. Interestingly, species richness has no direct impact on stability but does reduce asynchrony, indicating indirect complexity in the stability-dynamics relationship. This highlights the role of environmental and community factors in ecosystem stability, with asynchrony serving as a crucial intermediary in maintaining stability in fluctuating conditions.
Structural equation model (SEM) of the relationship between fundamental balance, asynchrony, richness, nutrients, temperature and temporal stability. The model shows that fundamental balance has a negative effect on asynchrony, which in turn has a positive effect on temporal stability. The model also shows that fundamental balance has a direct negative effect on temporal stability. Temperature has a direct negative effect on temporal stability, while nutrients have a direct positive effect on temporal stability. Richness has a direct negative effect on asynchrony, but no direct effect on temporal stability.
